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Abstract 

Metropolis Monte Carlo simulation is used to investigate the elasticity 
of torsionally stressed double-stranded DNA, in which twist and supercoiling 
are incorporated as a natural result of base-stacking interaction and backbone 
bending constrained by hydrogen bonds formed between DNA complementary 
nucleotide bases. Three evident regimes are found in extension versus torsion 
and/or force versus extension plots: a low- force regime in which over- and un- 
derwound molecules behave similarly under stretching; an intermediate-force 
regime in which chirality appears for negatively and positively supercoiled 
DNA and extension of underwound molecule is insensitive to the supercoil- 
ing degree of the polymer; and a large-force regime in which plectonemic 
DNA is fully converted to extended DNA and supercoiled DNA behaves quite 
like a torsionless molecule. The striking coincidence between theoretic calcu- 
lations and recent experimental measurement of torsionally stretched DNA 
[Strick et al., Science 271, 1835 (1996), Biophys. J. 74, 2016 (1998)] strongly 
suggests that the interplay between base-stacking interaction and permanent 
hydrogen-bond constraint takes an important role in understanding the novel 
properties of elasticity of supercoiled DNA polymer. 

Introduction 





Recent years have witnessed a remarkably intense experimental and theoretical activity in 
searching for the elasticity of a single supercoiled DNA molecule (see, e.g. Strick et al., 1996, 
1998; Fain et al., 1997; Vologodskii and Marko, 1997; Moroz and Nelson, 1997; Bouchiat 
and Mezard, 1998, Zhou et al., 1999). Within a cell, native double-stranded DNA (dsDNA) 
often exists as a twisted, and heavily coiled, closed circle. Differing amount of supercoiling, 
in addition to affecting the packing of DNA within cells, influences the activities of proteins 
that participate in processes — such as DNA replication and transcription — that require 
the untwisting of dsDNA (Wu et al., 1988). It is believed that changes in supercoihng can 
also promote changes in DNA secondary structure that influences the binding of proteins 
and other hgands (Morse and Simpson, 1988). 

In recent experiments (Strick et al., 1996, 1998) on single torsionally constrained DNA 
molecule, it was found that the supercoiling remarkably influences the mechanical property 
of DNA molecules. When apphed with relatively low stretching force, a supercoiled molecule 
can reduce its torque by writhing, forming structures known as plectonemes. Therefore, the 
distance between two ends of the polymer decreases with increasing supercoiling. But above 
a certain critical force fc, this dependence of extension on supercoiling disappears. More 
strikingly, the value of fc is significantly different for positively and negatively supercoiled 
DNA, i.e. fc ~ 0.8pN for underwound molecule and fc ~ 4.5pN for overwound ones. On 
the theoretical side, harmonic twist elasticity and bending energy according to the wormlike 
chain model have been used to understand supercoihng of DNA polymer (Fain et al., 1997; 
Vologodskii and Marko, 1997; Moroz and Nelson, 1997; Bouchiat and Mezard, 1998), and 
some qualitative mechanic features of plectonemic structures of supercoiled DNA polymer 
have been described by the harmonic twist elasticity. But because of the chiral symmetry 
of harmonic twist elasticity, the asymmetry of elastic behaviors of supercoiled DNA can not 
be understood by this model, and especially the three obvious mechanic regimes observed in 
experiment of supercoiling DNA (Strick et al., 1996, 1998) still need better understanding. 
To understand the supercoiling property as well as the highly extensibility of DNA, we 
have developed a more realistic model in which the double-stranded nature of DNA is taken 
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into account explicitly (Zhou et al., 1999). The supercoiling property of highly extended 
DNA was investigated analytically. Here, we aim at performing a thorough and systematic 
investigation into the property of supercoiled DNA by using Monte Carlo simulations based 
on this model. 

As we have known, the bending energy of DNA polymer is mainly associated with the 
covalent bonding between neighboring atoms of DNA backbone (Nossel and Lecar, 1991). In 
our previous work (Zhou et al., 1999), van der Waals interactions between adjacent basepairs 
was introduced and this helps to explain the highly cooperative extensibility of overstretched 
DNA (Cluzel et al., 1996; Smith et al., 1996). And it has been shown that the short-range 
base-stacking interaction takes a significant role in determining the elastical property of 
DNA. Lennard- Jones type potential between adjacent basepairs can be written as 



(see also Fig. 1). The folding angle 9 of the sugar- phosphate backbones around DNA central 
axis is associated with the steric distance r of adjacent basepairs by r = rocos6', where Vq 
is the backbone arclength between adjacent bases. The asymmetric potential related to 
positive and negative folding angle 9 in Fig. 1 ensures a native DNA to take a right-handed 
double- helix configuration with its equilibrium folding angle 6'eq ~ 6*0. This double-helix 
structure is anticipated to be very stable since e (~ 14 ksT) is much higher than thermal 
energy ksT according to the results of quantum chemical calculations (Saenger, 1984). 

In case that DNA polymer is torsionally constrained, the basepair folding angle will 
deviate from the equilibrium position ^eq- However, if the stretching force is very small, the 
folding angle may deviate from ^eq only slightly. This is because of the following reason: 
As we can infer from Fig. 1, the base-stacking potential is very sharp around 9q, and a 
relatively large force is needed to make 9 deviate considerably from its equilibrium value. 
It is reasonable for us to anticipate that a supercoiled DNA under low stretching force will 
convert its excess or deficit linking number into positive or negative writhing of its central 
axis. Since the central axis is symmetric with respect to positive or negative writhing, the 
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elastic response of DNA at this force regime will certainly be symmetric with positive or 
negative degree of supercoiling. Only when the stretching force becomes large enough will 
the chirality of supercoiled DNA appear. In this regime, it becomes more and more difficult 
for the central axis to writhe to absorb linking number and an increasing portion of the 
linking number will be converted to twisting number of the backbones, which will certainly 
changes the twisting manner of dsDNA. Since Eq. 1 shows that for dsDNA untwisting is 
much easier than overtwisting, chiral behavior is anticipated to emerge. This opinion is 
consistent with the experimental result of Strick et al. (1996). 

In this paper, we investigate the mechanical properties of supercoiled DNA by numerical 
Monte Carlo method. Base-stacking van der Waals interactions between adjacent basepairs 
are incorporated by introducing the new degree of freedom, namely the folding angle 6. 
A fundamental difference from the previous approaches (See, e.g., Vologodskii and Marko, 
1997), which try to include the twist degrees of freedom by adding extra terms to the free 
energy, is that the twist and supercoiling are treated as the cooperative result of base-stacking 
and backbones bending constrainted from permanent basepairs. The striking coincidence 
between theoretic calculations and experimental data of supercoihng DNA (Stick et al., 
1996, 1998) indeed confirms this treatment. 

Model and method of calculation 

In the simulation, the double-stranded DNA molecule is modeled as a chain of discrete 
cylinders, or two discrete wormlike chains constrained by basepairs of fixed length 2R (Fig. 
2). The conformation of DNA molecule of N straight cylinder segments is specified by the 
space positions of vertices of its central axis, rj = (a;(i), y(i), z{i)) in 3-D Cartesian coordinate 
system, and the folding angle of the sugar-phosphate backbones around the central axis, 
6i, i = 1,2, ■ ■ ■ , N. Each segment is assigned the same amount of basepairs, ntp, so that the 
length of the ith segment satisfies 

A I I . COS^j 

Asi = - rj_i = 0.34nbp7 (2) 

(cos6')o 
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where (• • •)o means the thermal average for a relaxed DNA molecule. Moreover, bearing in 
mind the experimental fact that there are about 10.5 basepairs for each turn of a native 
double helix DNA and the average distance between the adjacent basepairs is about do — 
0.34nm, we have set the basepair length as 2R — (10.5do/7i')(tan^)o in our model. 

Metropolis Monte Carlo method (Metropohs et al., 1953) is used to simulate the equi- 
librium evolution procedure of torsionally stretched dsDNA molecule. At each step of the 
simulation procedure, a trial conformation of the chain is generated by a movement from the 
previous one. The starting configuration is chosen arbitrarily (except that some topology 
and bound conditions should be satisfied, see below) and the averaged results of equilibrium 
ensemble are independent of the initial choice after numerous movements. The probability 
of acceptance of the movement depends on the difference in energy between the trial and 
the current conformations, according to the Boltzmann weight. When a trial movement is 
rejected, the current conformation should be counted once more. This procedure is repeated 
numerous times to obtain an ensemble of conformations that, in principle, is representative 
of the equilibrium distribution of DNA conformation. 

The DNA model 

As we have known, double strand DNA is formed by winding two polynucleotide back- 
bones right-handedly around a common central axis. Between the backbones nucleotide 
basepairs are formed with the formation of hydrogen bonds between complementary bases. 
In our continuous model (Zhou et al., 1999), the embeddings of two backbones are defined 
by ri(s) and r2{s'). The ribbon structure of DNA is enforced by having r2(s') separated 
from ri(s) by a distance 2R, i.e. r2(s') = ri(s) + 2i?b(s) where the hydrogen-bond-director 
unit vector b(s) points from ri(s) to ^2(3'). As the result of the wormhke backbones, the 
bending energy of two backbones can be written as 



The formation of basepairs leads to rigid constraints between the two backbones and at 
the same time they hinder considerably the bending freedom of DNA central axis because 
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of the strong steric effect. In the assumption of permanent hydrogen bonds (Everaers et 
al., 1995; Liverpool et al., 1998; Zhou et al, 1999), |s' — s| =0. The relative sliding 
of backbones is prohibited and the basepair orientation lies perpendicular to the tangent 
vectors ti = dri/ds and t2 = dr2/ds of the two backbones and that of the central axis, t: 
b ■ tl = b ■ t2 = b ■ t = 0. By defining the folding angle as half of the rotation angle from 
t2(s) to ti(s), i.e., the intersection angle between tangent vector of backbones ti(2) and DNA 
central axis t, we have 

{tl = cos 61 1 + sin 6* b X t 
(4) 
t2 = COS 6^ t — sin b X t. 

Therefore, the bending energy of the two backbones can be rewritten as 

where ds denotes arc-length element of the backbones, L the total contour length of each 
backbone, and n is the persistence length of one DNA backbone. Bearing in mind that 
the pairing and stacking enthalpy of the bases significantly increase bending stiffness of 
polymer axis, the experimental value of persistent length of dsDNA polymer is considerably 
larger than that of a DNA single strand (See, e.g. Smith et al., 1996). To incorporate 
the steric effect and also considering the typical experiment value of persistent length of 
dsDNA p = 53nm, the simpliest way is to substitute k in the first term of Eq. |^ with a 
phenomenological parameter k* = 53.0/ 2 {cos 9) oiamkBT (Zhou et al., 1999), hereafter this 
is assumed. 

Taking into account Eqs. |l| and |^, the total energy of dsDNA molecule with segments 
in our discrete computational model is expressed as 

N-l TV -^bp 

^ = a E 7f + a' E (^^+1 - + m E sin3 9, tan 9, + J2 U{9,) - fz{N), (6) 

i=l 1=1 ^ i=l j=l 

where 7^ is the bending angle between the {i — l)th and the ith segments (Fig. 2), Abp the 
total number of basepairs of DNA polymer, and z{N) is the total extension of the DNA 
central axis along the direction of the external force / (assumed in the 2;-direction). 
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Since Kuhn statistical length of dsDNA polymer is associated with its bending stiffness 
(the Kuhn length is twice as persistence length of dsDNA polymer according to the wormlike 
chain model), one can decide bending rigidity parameter a of the discrete chain accordingly. 
Suppose that we take the N discrete segments to simulate the behaviors of a dsDNA polymer 
of n Kuhn statistical length, the length of m(= N jn) segments should correspond to one 
Kuhn statistical length. Therefore, for any chosen value m, we can decide the bending 
rigidity parameter a in the way (see Appendix) 

1 + (cos 7) 



1 — (cos 7) ' 
where 



(7) 



/ \ /(7cos7exp(-a72)sin7ci7 

cos 7 = — — . (8) 

Jq exp(— 0:7^) sm 7(27 

In principle, the discrete DNA model becomes continuous only when m approaches 
infinity. The CPU time needed for a simulation, however, increases approximately as 
A^^ = innif' . So it is necessary to choose a value of m that is large enough to ensure 
reliable results but small enough to keep the computational time within reasonable bounds. 
Our calculation and also previous work (Vologoskii et al., 1992) showed that simulated 
properties do not depend on m if it exceeds 8. Therefore, m = 8 was used in the current 
calculation, for which the bending constant a = l.SdbksT. Furthermore, we have chosen 
A^ = 160 in consideration of the feasible computer time. Since Kuhn statistical length of 
dsDNA is taken as 106nm, the B-form length of the polymer in our simulation corresponds 
to Lb = 2120nm or 6234 base-pairs. The constant a' in the second term of Eq. ^ should 
be associated with stiffness of the DNA single strand. As an crude approximation, we have 
taken here a' = a = 1.895fcBT.Q 

The fourth term in Eq. |^ accounts for van der Waals interactions between adjacent 
basepairs (see Eq. |l]). 9o {= 62°) is related to the equilibrium distance between a DNA 



^ Our unpublished data show that, the amount of second term of Eq. |^ is quite small compared 
with other four terms. And the result of simulation is not sensitive to a' . 
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dimer. The base-stacking intensity e is generally influenced by composition and sequence 
of nucleotide chains. For example, the solubility experiments in biphasic systems show that 
stacking interactions between purine and pyrimidine bases follow the trend 

purine — purine > pyrimidine — purine > pyrimidine — pyrimidine. 

Since we do not distinguish the specific base-sequence of purine and pyrimidine in our DNA 
model, we take statistic average of stack energies as e = 14A;bT, according to the result of 
quantum chemical calculations (Saenger, 1984). 

To simulate the extension of the stretched DNA chain, we fixed one of its ends at original 
point in 3-D Cartesian system and applied a force / directed along the z axis to the second 
end, which corresponds to the fifth term of Eq. |^. 

Calculation of link number 

The number of times the two strands of DNA double helix are interwound, i.e., the link 
number Lfc, is a topologic invariant quantity for closed DNA molecule and also for linear 
DNA polymer in case that the orientations of two extremities of the linear polymer are 
fixed and any part of polymer is forbidden to go round the extremities of the polymer. An 
unstressed B-DNA molecule has one right-handed twist per 3.4nm along its length, i.e., 
Lfco = Lb/3.4. Under some twist stress, the link number of DNA polymer may be different 
from its torsionally relexed value. In all case when AL/c = Lk — Lk^ ^ 0, the DNA polymer 
is called "supercoiled" (Vologodskii and Cozzarelli, 1994). The relative difference in link 
number 

Lk - Lko 

signifies the degree of supercoiling which is independent upon the length of DNA polymer. 
The native DNA of organisms living at physiological environment are found always slightly 
underwound and its supercoiling degree is between —0.03 and —0.09 (Bauer, 1978; Volo- 
godskii and Cozzarelli, 1994), which is believed significantly relevant in some fundamental 
biological processes (Wu et al., 1988; Morse and Simpson, 1988). 
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In addition to counting directly the number of times the two strands are interwound, 
the hnk number of closed DNA circle can be conveniently calculated by White's theorem 
(White 1969) 



The twist Tw is the number of times basepair twist around central axis and does not depend 
upon the configuration of molecule axis. The writhe Wr of molecule is a simple function of 
only the molecule axis vector r(s) (White, 1969; Fuller, 1971) 



Wr is scale invariant and dimensionless and changes sign under reflection or inversion of r, 
reflecting the cross product in the formula above. Therefore Wr = if r(s) is planar or 
otherwise reflection symmetric. 

In order to control and measure experimentally the supercoiling degree of linear DNA 
polymer, Strick et al. (1996, 1998) attached one end of DNA molecule to a glass cover slip 
by DIG-anti-DIG links and other end to a paramagnetic bead by biotin-streptavidin links. 
Bearing in mind the diameter of magnetic bead (~ 4.5/im) is far beyond that of polymer, 
the anchoring points can be considered as on impenetrable walls and ~ 16-/im-long DNA 
(Strick et al., 1996) in fact is prohibited to pass around the ends of the polymer. A magnetic 
field pointing in the plane of the microscope slide was applied to fix the orientation of the 
bead. Therefore, by rotating the magnets and counting the time of turns, the link number 
Lk of the linear DNA molecule can be controlled and measured experimentally. 

In Monte Carlo calculations, we restrict the DNA chain by two impenetrable parallel walls 
crossing the chain ends which is to simulate the above mentioned experimental equipment 
of the magnet bead and the microscope slide (see also the treatment in Vologodskii and 
Marko, 1997). The walls are always parallel to xy plane in our Cartian coordinate system 
and thus perpendicular to the direction of the force applied to the chain ends. 

One way to calculate the link number Lk of DNA molecule in our Monte Carlo simulation 
is to use the White's formula Eq. [10|. However, the writhe Wr is defined only for closed 



Lk = Tw + Wr. 



(10) 




(11) 
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chain. In order to solve the problem, we add three long flat ribbons to the two ends of 
the DNA chain in each conformation during the simulation procedure. The axes of these 
ribbons are kept in the same planar and consist a closed circle together with the linear DNA 
chain. Since there is no any twist in the added three flat ribbons, it is not difficult to verify 
from Fig. 3 that the number of times two strands interwind Lki in Fig. 3a is equal to the 
link number of new closed polymer Lkc in Fig. 3d. Therefore, we only calculate Lk of the 
closed chain in our simulations according to Eqs. ^ and |ll]. 

Quite similar to the model by Tan and Harvey (1989) in which the twist of each base-pair 
of DNA chain is explicitly specified, the folding angle of backbones in each segments has 
been given in our model. So the twist can be directly calculated by 

1 ^ 

Tw = — -yAsitanOi. (12) 
The writhe Wr of the new DNA circle can be calculated according to Eq. |Tl|. 
Simulation procedure 

For any given force, equilibrium sets of conformations of DNA chain are constructed 
using the Metropolis MC procedure (Metropolis et al., 1953). Three kinds of movements 
have been considered in our simulations (see Fig. 4). 

In the first type of movement, a random chosen segment is undertwisted or overtwisted 
by an angle Ai. In other words, the folding angle 6i of the chosen segment is modified into 
a new value 9'^ = 9i + Xi- When 9'- is beyond the setting interval [—9^,9^] from one side, 
it will re-enter the interval from the opposite side according to the periodicity assumption. 
Although the geometric limit of folding angle of DNA backbone is 9^ = v/2, we set 9^ = 85° 
here to avoid the possible divergency in numerical calculation of potential of Eq. |I|. It should 
be mentioned that, this movement modifies not only the folding angle of the chosen segment 
but also the coordinates of all the behind vertices rj,j = i,---,N along the length, since 
when the folding angle 9i is changed we have also changed the length of the segment Asi 
according to Eq. ^. So we should translate all those segments behind this one to make the 
chain match up (Fig. 4a). 
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In the second type of movement, an interval subchain containing arbitrary amount of 
segments will be rotated by an angle of A2 around the straight line connecting the vertices 
bounding the subchain (Fig. 4b). The third type of movement involves a rotation of the 
subchain between a chosen vertices and the free end by an angle of A3, around an axis with 
arbitrary orientation (Fig. 4c). All three types of movements satisfy the basic require- 
ment of the Metropolis procedure of microscopic reversibility, i.e. the probability of trial 
conformation B when current conformation is A must be equal to the probability of trial 
conformation A when current conformation is B. 

All three types of movements change the configurations of DNA chain. But from the 
viewpoint of energy, their functions are quite different. While the first type of movement 
concerns mainly with modifying twist and stacking energy, the second one changes only the 
bending energy and the third modifies both bending energy and extension of DNA chain. 
Each of them is performed in the probability of 1/3. The value of Ai, A2, A3 are uniformly 
distributed over interval (-A?,A?), (-A^,A^) and (-A^, Ag) respectively, and A?, A^ and A^ 
are chosen to guarantee that about half of the trial moves of each type are accepted. 

The starting conformation of DNA chain is unknotted. But the configurations after 
numerous steps of movements may become knotted, which violates the topologic invariance 
of chain and is incorporeal. Especially, both ends of molecule are anchored in the experiment 
and knots never occur. In order to avoid knotted configuration, we should check the knot 
status for each trial conformation. The most effective way to clarify the knot categories of 
DNA circle is to calculate its Jones polynomial (Jones, 1985), which is strictly topological 
invariant for knot categories. But the computational calculation of Jones polynomial is quite 
prolix at this moment. In our case that it is only necessary to distinguish between unknot 
and knot categories, the classical Alexander polynomial (Alexander, 1928; Conway, 1969) 
is enough to meet this requirement although it is of weaker topological invariants and does 
not distinguish mirror images. For trivial knot, Alexander polynomial A{t) — 1; and A{t) is 
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usually not equal to 1 for knotted chain. [| Convenient algorithms for computer calculation 
of Alexander polynomial had been well built (see, e.g. Vologodskii et al. 1974; Harris and 
Harvey, 1999). We only calculate the value of A(— 1) in our simulation. In case that the 
trial movement knots the chain, the energy of trial conformation is set to be infinite, i.e. it 
will be rejected. 

Another interaction considered in our simulation is the steric effect of polymer chain. 
Since the segment has finite volume, other segments cannot come into its own space region. 
This interaction evidently swells the polymer (Doi and Edwards, 1986). To incorporate 
this exclude-volume effect into our simulation, for each trial conformation, we calculate the 
distance of between any point on the axis of a segment and any point on the axis of another 
non- adjacent segment and check whether this distance is less than the DNA diameter 2R. 
If the minimum distance for any two chosen segments is less than 2R, the energy of trial 
conformation is set infinite and the movement is rejected. 

During the evolution of DNA chain, the supercoiling degree a may distribute around 
all the possible values. In order to avoid the waste of computation events, we bound the 
supercoiling a of DNA chain inside the experimental region (Strick et al., 1996, 1998), i.e. 
—0.12 > a > 0.12. When the torsion degree of trial conformation is beyond the chosen 
range, we simply neglect the movement and reproduce a new trial movement again. 

Result of Monte Carlo simulation 

To obtain equilibrium ensemble of DNA evolution, lO'' elementary displacements are 
produced for each chosen applied force /. The relative extension x and supercoiling degree 
cr of each accepted conformation of DNA chain are calculated. When the trial movement is 
rejected, the current conformation is count up twice (see Metropolis et al., 1953). 

^ Although there are nontrivial knots whose Alexander polynomials equal unity, this case is very 
rare. One of the example for nontrivial knot with A(t) = 1 can be found in Vologodskii et al. 
(1974). 
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In order to see the dependence of mechanics property of DNA upon supercoihng degree, 
the whole sample is partitioned into 15 subsamples according to the value of the supercoihng 
degree a. For each subsample, we calculate the averaged extension 

(^,) = 1^E£i^, J = 1,---,15 (13) 

and the averaged torsion 

M = ^E^- i = l,---,15, (14) 

^^j i=i 

where Nj is the number of movements supercoihng of which belong to jth subsample. 

We display the force versus relation extension for all positive and negative supercoihng 
in Fig. 5a and c respectively. As a comparison, the experimental data (Strick et al., 1998) 
are shown in Fig. 5b and d. In Fig. 6 is shown the averaged extension as a function 
of supercoihng degree for 3 typical applied forces. At low force, the extension in our MC 
simulation saturates at a value greater than zero because of the impenetrable walls which 
astrict the vertical coordinate of the free end always higher than that of any other points of 
the DNA chain. The same effect of the impenetrable walls was found in earlier works (see 
Fig. 9 of the paper by Vologodskii and Marko, 1997). For conciseness, we did not show the 
points the relative extension of which is less than 0.15 in Fig. 5 and 6. 

In spite of quantitative difference between Monte Carlo results and experimental data, 
the qualitative coincidence is striking. Especially, three evident regimes exist in both exper- 
imental data and our Monte Carlo simulations: 

i). At a low force, the elastic behaviour of DNA is symmetrical under positive or negative 
supercoihng. This is understandable, since the DNA torsion is the cooperative result 
of hydrogen-bond constrained bending of DNA backbones and the base-stacking in- 
teraction in our model. At very low force, the contribution from applied force and 
the thermodynamic fluctuation perturbate the folding angle 9 of basepair to derive 
very little from the equilibrium position ^o- Therefore, the DNA elasticity is achiral 
at this region (see the Introduction part of this paper). For a fixed applied force, the 
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increasing torsion stress tends to produce plectonemic state which shorten the distance 
of two ends, therefore, the relative extension of hnear DNA polymer. These features 
can be also understood by the traditional approaches with harmonic twist and bending 
elasticity (Vologodskii and Marko, 1997; Bouchiat and Mezard, 1998). 

ii) . At intermediate force, the folding angle of basepairs are pulled slightly further away 

from equilibrium value Oq where van der Waals potential is not symmetric around 
Oq. So the chiral nature of elasticity of the DNA molecule appears. In negative 
supercoiling region, i.e. 9 < 9q, the contribution of applied force dominates that 
of potential because of the low plateaus of U{9). So the extension is insensitive to 
negative supercoiling degree. On the other hand, the positive supercoiling still tends 
to contract the molecule. 

iii) . At higher force, the contribution of the applied force to the energy dominates that of 

van der Waals potential in both over- and underwound DNA. The extension of DNA 
accesses to its B-form length. Therefore, the plectonemic DNA is fully converted to 
extended DNA, the writhe is essentially entirely converted to twist and the force- 
extension behaviour reverts to that of untwisted (cr = 0) DNA as expected from a 
torsionless worm-hke chain model (Smith et al., 1992; Marko and Siggia, 1995; Zhou 
et al., 1999). Because of the effect of impenetrable wall, however, the extension of 
DNA molecule in our calculation is slightly higher than experimental data. 

In conclusion, the elasticity of supercoiled double-stranded DNA is investigated by Monte 
Carlo simulations. In stead of introducing an extra twist energy term, twist and supercoiling 
are leaded into as a nature result of cooperative interplay between base-stacking interac- 
tion and sugar-phosphate backbones bending constrained by permanent hydrogen bonds. 
Without any adjustable parameter, the theoretic results on the correlations among DNA 
extension, supercoiling degree and applied force agree qualitatively to recent experimental 
data by Strick et al (1996, 1998). 

It should be mentioned that there is an up-limit of supercoiling degree for extended 
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DNA in current model, i.e. CTmax ~ 0.14, which corresponds to ^ = 90" of folding angle. 
In recent experiments, AUemand et al. (1998) twisted the plasmid up to the range of 
— 5 < cr < 3. They found that at this "unrealistically high" supercoiling, the curves of force 
versus extension for different a split again at higher stretch force (> 3pN). As argued by 
AUemand et al. (1998), in the extremely under- and overwound torsion stress, two new 
DNA forms, denatured-DNA and P-DNA with exposed bases, will appear. In fact, if the 
deviation of the angle which specifies DNA twist from its equilibrium value exceeds some 
threshold, the corresponding torsional stress causes local distraction of the regular double 
helix structure (Vologodskii and Cozzarelli, 1994). So the emergence of these two striking 
forms is essentially associated with the broken processes of some basepairs under super- 
highly torsional stress. In this case, the permanent hydrogen constrain will be violated 
and the configuration of base stacking interactions be varied considerably. We hope, with 
incorporation of these effects at high supercoiling degree, our model should reproduce the 
novel elastic behaviour of DNA. This part of work is in progress. 
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Appendix: Kuhn Statistical Length of Discrete Chain 

Let us consider a discrete chain of segments with each of length /q, the end-to-end 
vector of which is written as 



N 



R = ^0 



(15) 



i=l 



where tj = 



|Ri-Ri_i|- 
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For chains with bending stiffness, e.g. the DNA model described in Eq. (tj+fc ■ tj) does 
not vanish for /c 7^ 0. tj+fc can be expressed relatively to i + /c — I'th segment as 

ti+fc = cos7i+fc_itj+fe_i + sin7i+fc_ini+fc_i, (16) 

where ■ji+k-i is the bending angle between i + A; — I'th and i + fc'th segments as defined 
in Eq. and nj_|_fc_i is the unit vector coplanar with tj+fc and tj+fc_i but perpendicular to 
the latter. If the average of tj+fc is taken with the rest of the chain (i.e., tj, tj+i, • • • , tj+fc_i) 
fixed, one obtains 

(tj+A,.)ti,ti+i,---,t,+fe_ifixcd = (cOS7i+fe_i)ti+fc_i, (17) 

since {ni+k-i)u,u+i,-,ti+k-ifi^ed = according to Eq.|. Multiplying both sides of Eq. |T3 by 
tj and taking the average over tj, t^+i, • ■ • , tj+fc_i, one has 

(tj+fc ■ ti) = (cos7)(tj+fc_i ■ ti), (18) 

where (cos 7) is not specific to segments and given by Eq. |^. This recursion equation, with 
the initial condition = 1, is solved by 

(ti+fc-t,) = (cos7)^ (19) 

Thus for large N, (R^) is given by 

TV TV 

i=i j=i 

TV-l N-i 

= /g(iV + 2^ ^(t,-t,+,)) 

i+l k=l 

2I + (cos 7) 



1 — (cos 7) 

Therefore, Kuhn statistical length of the discrete chain can be written as 

„^<m.,,i±pA (20, 

-Rmax 1-(C0S7) 

where -Rmax is the maximum length of the end-to-end vector. 
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FIGURES 

FIG. 1. The van der Waals interaction potential versus folding angle of sugar-phosphate 

backbones around DNA molecule axis. 

FIG. 2. The configuration of discrete DNA chain in our model. 

FIG. 3. The schematic diagram to calculate link number in our simulations, (a). For a linear 
supercoiled DNA chain with one end attached to a microscope slide and with another end attached 
to a magnetic bead, when the orientation of the bead is fixed and the DNA chain is forbidden to go 
round the bead, the number of times for two strands to interwind each other, the linking number of 
the linear DNA (Lki), is a topological constant, (b). The DNA double helix is stretched to a fully 
extended form while the orientation of bead keeps unchanged. The link number of linear DNA 
chain is equal to the twist number, i.e. Lki = Twi. (c). Three long flat ribbons are connected to the 
two ends of the linear twisted DNA of (b) . The link number of the new double helix circle is equal 
to that of linear DNA chain, i.e. Lkc = Twc = Twi = Lki since the writhe of the rectangle loop is 
0. (d). The DNA circle in (c) can be deformed into a new circle, one part of which has the same 
steric structure as the linear supercoiled DNA chain in (a). So by adding three straight ribbons, 
the link number of linear double helix DNA can be obtained by calculating the link number of the 
new DNA circle, i.e. Lki = Lkc = Tw + Wr. 

FIG. 4. Trial motions of the DNA chain during Monte Carlo simulations. The current confor- 
mation of DNA central axis is shown by solid lines and the trial conformation by dashed lines, (a). 
The folding angle in ith segment Oi is changed into 9i + Xi. All segments between ith vertex and 
the free end are translated by the distance of |Asj — As^|. (b). A portion of the chain is rotated by 
an angle of A2 around the axis connecting the two ends of rotated chain, (c). The segments from a 
randomly chosen vertex to the free end are rotated by an angle A3 around an arbitrary orientation 
axis which passes the chosen vertex. 
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FIG. 5. Force versus relative extension curves for negatively (a,b) and positively (c,d) super- 
coiling DNA molecule. Left two plots (a) and (c) are the results of our Monte Carlo simulation, 
and the horizontal bars of points denote the statistic error of relative extension in our simulations. 
Right two plots (b) and (d) are the experimental data (Strick et al., 1998). The solid curves serve 
as guides for the eye. 

FIG. 6. Relative extension versus supercoiling degree of DNA polymer for three typical stretch 
forces. Open points denote the experimental data (Strick et al., 1998) and solid points the results 
of our Monte Carlo simulation. The vertical bars of the solid points signify the statistic error of 
the simulations, and the horizontal ones denote the bin-width that we partition the phase space of 
supercoiling degree. The solid lines connect the solid points to guide the eye. 
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